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The matter created in relativistic heavy ion collisions is fairly well described by ideal 
hydrodynamics, and somewhat better described by viscous hydrodynamics. To this point, 
most viscous calculations have been two-dimensional, based on an assumption of Bjorken 
boost invariance along the beam axis. Here, first results are presented for a fully three- 
dimensional viscous model. The model is described and tests of the numerical accuracy 
of the code are presented. Two- and three-dimensional runs are compared, and modest 
changes are observed for mid-rapidity observables at the highest RHIC (Relativistic Heavy 
Ion Collider) energies. 



I. INTRODUCTION 



Hydrodynamics has been central to the interpretation of experimental results from the Relativistic 
Heavy Ion Collider (RHIC) [HE]- Even relatively simple hydrodynamic models that include neither 
the effects of viscosity, a realistic phase transition, nor hadronic rescattering are able to roughly 
reproduce such important observables as the transverse momentum spectra and the large anisotropic 
flow observed at mid-rapidity [3]. Despite difficulties in describing these mid-rapidity observables 
simultaneously with the femtoscopic source sizes and elliptic flow at non-zero longitudinal rapidity 
[H [5] , the results strongly suggested underlying collective behavior and were taken as strong evidence 
for the formation of a Quark-Gluon Plasma (QGP). 

Since that time, considerable effort has been expended by the modeling community to characterize 
and constrain the properties of the matter created in these heavy ion collisions. In particular, several 
models have focused on developing viscous hydrodynamics for the detailed description of anisotropic 
flow at mid-rapidity to constrain the viscosity of the QGP [6, 7J. To simplify the calculations, most 
of these models have restricted themselves to mid-rapidity and used the expected approximate boost 
invariance near zero rapidity [8] to motivate two-dimensional treatments. 

In reality, the spectrum of produced particles is well described by a Gaussian of radius wl.6 units in 
rapidity 0, and one questions the extent to which the finiteness of the extent affects the evolution at 
mid-rapidity. Furthermore, the character of the results away from mid-rapidity display some interesting 
trends that likely provide some unique and useful constraints on the character of the matter produced 
and the important features of the theory that governs the production of that matter. To this end, we 
endeavor to explore the predictions of viscous hydrodynamics in the full configuration space. 



II. HYDRODYNAMIC THEORY 

We consider a system of interacting particles. When the density of particles is relatively low and 
the mean free paths for particles in the systems are larger than characteristic inter-particle separations, 
the system can be described as a series of binary collisions. One can then describe the evolution with a 
Boltzmann equation, which requires following the evolution of the six-dimensional phase space density, 
/(p, r, t), for each species. When the mean free path is much smaller than the characteristic size of the 
system, the local phase space distribution takes on a thermal form, and the system and its evolution 
can be described by following a small number of variables in configuration space, i.e., current densities 
(p^^) and the stress-energy tensor (T^u)- The equations of motion for the quantities are driven by the 
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conservation laws: 



d^.P^[x) = 0, (1) 
d^T^^^x) = 0. (2) 

In terms of thermodynamic variables, i.e. energy density (e) and pressure the stress-energy 
tensor for ideal hydrodynamics is 

T^"' = (e + - g^'^'p (3) 

where g^^'^ = diag(l, -1, -1, -1) is the metric tensor. This expression makes clear that one still requires 
an equation of state to link the thermodynamic variables - i.e. p(e, p) - in order to close the system of 
equations. 

If the distribution function is allowed to deviate modestly from the equilibrium value, the theory 
may be extended to account for the effects of a finite mean free path [3 [TU] with viscous treatments. 
Roughly, one assumes that collisions cause the stress-energy tensor to relax toward the Navier-Stokes 
value or that the corresponding modification to the stress energy tensor (vr^'^) leads to a quadratic 
increase to the equilibrium entropy density [Tl] . The equations of motion are then taken either from 
higher moments of the Boltzmann equations or from the local conservation of the modified entropy 
density jl2H16j . Such a construction requires the introduction of several new transport coefficients, 
but supplies equations of motion that apply to systems further from equilibrium. In this case, we will 
only allow the transport coefficients associated with shear viscosity to be non-zero. This introduces 
two new quantities: the shear viscosity (r/) and the shear relaxation time (t^). The resulting equations 
of motion are the conservation of stress energy tensor and the relaxation of the shear tensor to the 
Navier-Stokes value [T2HIi] : 

T^'^ = (e + p)u^'u^ - g^^p + n"" (4) 

= -rj (d'n^ + B^n' - (2/3) (9 • n)6i,^ . (6) 

The tildes denote quantities determined in the frame of the matter, and roman indices denote that 
only spatial indices are being considered. Here, vr^'^ is the viscous correction to the stress-energy 
tensor. In the frame of the matter 

^00 ^ 

= due its orthogonality to the collective velocity. The 
quantity cr^ is related to the equilibrium fluctuation of vTy, and has the same units of pressure [12]. 
For a conformal theory it can be replaced with the energy density |16j . 



For reasons of clarity, we choose to tabulate moments of the local shear tensor defined 

ai = \ - %yy) ■ a2 = (tt"" + - 2^^) ■ (7) 
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6 = (1/3) [t,^^ + T:yy + 

This choice of moments clearly separates the shear effects of the strong initial longitudinal expansion 
into the variable 02, whereas ai will contain only the effects of the anisotropic transverse expansion. 
Furthermore, we have limited the number of variables to the minimum - gaining speed but sacrificing 
the possibility of a check on the numerical accuracy of the solution. 

Hydrodynamic models best apply to the matter at small spatial rapidity soon after the rapidly 
receding nuclei have passed through one another. The spatial rapidity is defined by 

^s = \ln{^-^^; T = Vt^, (8) 

and is a measure of the longitudinal position, i.e., 6rjs = 5z/t. For a boost-invariant system, the 
collective velocity is Vz = z/t, which makes the spatial rapidity equal to the rapidity of the matter. 
The condition Vz = z/t describes a system where each density element had a trajectory that began at 
z = t = 0, and did not accelerate from that point. From the perspective of hydrodynamics, the lack of 
acceleration is expected if intrinsic quantities, such as energy density, vary slowly with rjs and depend 
mainly on the proper time r. This approximate symmetry motivates using hydrodynamic prescriptions 
where rjg is neglected, even in the presence of viscosity, which then reduces the dimensionality of the 
model from three to two. Of course, such a model only applies to observables measured near mid- 
rapidity. A principal goal of this study is to discern the validity of this approximation by comparing 
two- and three-dimensional versions. 

III. INITIAL CONDITIONS 

While considerable progress has been made into viscous hydrodynamics for event-by-event con- 
ditions with fluctuations |17H19j . for this paper we choose a simpler approach and assume averaged 
initial conditions such as those from the optical Glauber model [20]. We begin with a Woods-Saxon 
distribution of nucleons 

P^'^ = l+expKr-i?)/x] 

where R = 6.37 fm, x = 0-54 fm, and po = 0.17 fm^^ for Au nuclei. One then considers the total 
probability of interaction as the nuclei pass completely through one another, which is computed via 



the nuclear thickness function 

/oo 
dzp{x,y,z) (10) 
-oo 

The initial shape of the system is then scaled either to the number of wounded nucleons or the 
number of binary collisions. The difference between the scalings comes from the eligibility for multiple 
interactions. In the binary collision picture, nucleons are allowed to interact with one another multiple 
times during the crossing. This leads to scaling like the product of the two thickness functions 

UBcix, y, b) = aTA{x + 6/2, y)TB{x - 6/2, Y) (11) 

where a is the total inelastic cross-section for nucleon-nucleon pairs at this collision energy, Ta,b is 
the thickness function for a nucleus with A nucleons, and b is the impact parameter. 

In the wounded nucleon picture, pairs of interacting nucleons are removed via a combinatorial 

factor leading to 



nwNix, y, b) = Ta{x + 6/2, y) 



aTB{x-b/2,y) ^ ^' 
B 



+ sym. (12) 



where the symmetry is between the nuclei A and B. At lower (SPS) energies, it was found that the 
number of produced particles scaled with the number of wounded nucleons. For RHIC collisions, 
one expects that there may be deviation from this scaling due to the increasing importance of hard 
processes. Initial conditions are thus scaled partially to the number of wounded nucleons and partially 
to the number of binary collisions via 

e{x, y, b) = KQ [anwN{x, y, 6) + (1 - a)nBc{x, y, b)] , (13) 

where a is then the fraction of the initial condition scaled to the number of wounded nucleons. 

We then take this total number of interactions to be proportional to the initial energy density of the 
system at some time shortly after crossing. Both the overall energy scale (kq) and the thermalization 
time (ro) will be considered free parameters of the model, in addition to the fraction of the initial 
condition scaled to the number of wounded nucleons (a). In this paper we take tq = 0.8 fm/c, a = 0.85, 
and for the chosen value of kq, the central energy density at 6 = 2.21 fm is « 19.5 GeV/fm^. The 
choice was made to have the energy scale with n, rather than entropy, because energy density more 
justifiably scales with the number of independent superimposed sources. These parameters are chosen 
arbitrarily for the current model, as their values will be the subject of a future publication. 

The optical Glauber prescription sets the shape of the initial energy density in the transverse 
plane, it does not lend any insight into the longitudinal initial conditions. For that we turn to 
the experimental data, which shows little deviation from a simple Gaussian profile in longitudinal 



rapidity. Since Gaussians tend to remain roughly Gaussian in a hydrodynamic expansion |21] . we 
choose a Gaussian with a radius to be determined along with the other parameters of the model. For 
the results presented here, this radius was chosen to be 1.4 units of spatial rapidity. This choice of 
longitudinal source is in stark contrast with the long flat region with steep half-Gaussian used in other 
models [UESHSl]. 

Furthermore, while the Glauber model has proven useful in describing the initial state of the matter, 
a competing model (CGC) based on the assumption that the low momentum particles at midrapidity 
are generated by color fields has also been considered by many modelers |25H28j . Generally, CGC 
initial conditions have larger initial eccentricities and larger gradients in the tails of the transverse 
energy density, thus requiring somewhat larger shear viscosity to match experimental results [6]. 

With the energy density set, there remain nine components of the initial condition of the stress 
energy tensor to be set. In determining the initial condition for the collective velocity, hydrodynamic 
models have taken the approach that given a system which is newly thermalized, one might expect 
no transverse collective flow to be present [5l El [221124^ I29j . However, since some finite evolution of 
the system occurs after the crossing but before thermalization, it is unlikely that no evolution of the 
stress energy tensor should occur. If Bjorken boost invariance is valid, and if the dynamics are driven 
by a traceless stress energy tensor, one can predict the initial transverse acceleration regardless of the 
underlying theory [30]. For small r we assume the initial flow is a fraction ^ of this result. We will 
utilize the result 



where ^ is introduced as a free parameter of the model. In this paper we take ^ = 0.5, or half of the 
pre-equilibrium flow predicted. In the longitudinal direction, we will assume that the matter is not 
moving relative to the Bjorken expansion. 

Finally, the initial values of the shear tensor must be set. In the Landau frame, the shear tensor 
is constrained by its orthogonality to the fluid velocity (u^vr^'^ = 0), traceless-ness [g^^-K^^^ = 0), and 
symmetry (tt^'' = ■n'^^). This leaves five components to be determined. In recent history, other models 
have either set all the components to zero initially [22', W\:\ or set all the components to their Navier- 
Stokes values [7]. While using the Navier-Stokes values have been more prevalent, this approach leads 
to increasing values of [k^^ / P]{x) for large x, because rj/P cc at low temperature. Even for 
moderately small viscosities {rj/s ~ 0.16) and relevant temperatures (T ~ 100 MeV), one finds values 
of \tt^^\ that approach not just the pressure but the energy density. If e+p + n^^ = the conservation 
of longitudinal momentum approaches singularity and the longitudinal fluid velocity grows without 
bound. For initial conditions that reproduce experimental results, these problems tend to occur well 




(14) 
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outside of the freezeout surface. Also, since the relaxation times are very large when the viscous 
corrections are large, the evolution will tend not to approach such singularities if the initial conditions 
are not close to them. Because relaxation times tend to be short inside the freezeout surface, the 
system loses its memory of the initial values of the shear tensor within ~ 1.0 fm/c of the onset of 
hydrodynamics. For this study, the initial value of vr^^/P is chosen independent of position, leaving 
the effects of varying the parameter as a topic for future study. 



IV. ALGORITHM DESCRIPTION 

As discussed above, the model is designed to explore predictions of (3+l)-dimensional viscous 
hydrodynamics using averaged initial conditions. Since this avoids encountering the small-scale fluc- 
tuations of event-by-event hydrodynamics, the code does not invoke any aggressive scheme for solving 
the conservation equations in the presence of singularities or shocks. Instead, the ten equations of 
motion (Appx. B) are solved simultaneously using the second order Runge-Kutta (predictor-corrector) 
method on a fixed mesh in the expanding, Bjorken coordinates. Spatial derivatives of the integration 
variables are determined by second order central finite differences, and are corrected for the collective 
flow and the Bjorken expansion (Appx. A). 

The set of variables to be integrated are: the energy density, the fluid velocity relative to the 
Bjorken expansion, and the five components of the shear tensor in the matter frame, Oj. One of 
the advantages of choosing matter frame variables is that restrictions to the maximum deviation from 
equilibrium are more easily applied, though such restrictions are not used in this paper. The algorithm 
verifies the conservation the integral of T^^ in the mesh frame: 



= / (fxd^T^'^ = / drdrjdxdy 



,0. , roo + T- 



C = I T^^Tdrjdxdy + / dr 



T^'ds'+ I ^TdT]dxdy 



(15) 
(16) 



where the second term in Eq. 16 is the outgoing flux of T"* through a surface element ds*, the second 
term is the work done by the expansion of the mesh, and C is a constant. This quantity is conserved 
at the level of a tenth of percent in the three-dimensional viscous case for the cell densities used in 
this paper, which are dx = dy = dr]s = 0.15 (fm) and dr = 0.05 fm/c. A simple example of further 
investigation into the veracity of the code is provided in Appx. C. 



V. EQUATION OF STATE 



Significant progress on the transition region equation of state has been made since the early days 
of hydrodynamics |31H33j . Studies of lattice predictions have shown that adjusting for extrapolation 
to physical masses leads to a less dramatic phase transition at zero baryon chemical potential than 
originally assumed. This leads to an equation of state that is stiffer (higher speed of sound) through 
the transition region. The results presented use the three-flavor equation of state from the Wuppertal 
group [31j, but fits to MILC collaboration [33j data have also been explored. 

Of vital importance to the accuracy of transition to the hadronic cascade will be the matching of the 
energy density at the freeze-out temperature. While the above equation-of-state fits roughly reproduce 
data in the hadronic region, the energy density can disagree on the order of 10% with computations 
from the hadronic cascade itself. To ensure that this does not occur, we use the equation of state 
extracted directly from the hadron gas below the freeze-out temperature (T/), the lattice data above 
some temperature (Ti), and in between will merge the two. We merge continuously in the unitless 
entropy density ct(T) = s/T^ via the weighted linear function, a = wai + (1 — w)ah with weight 



«.(T) = i 



tanh I tan 



vr / Ti-T 

- 2— 1 

2 V Ti-Tf 



+ 1 



(17) 



which ensures the continuity and smoothness of cr(T) and therefore a continuous speed of sound 
squared via 



'^^ a dT 

We take T; = 190 GeV and Ty = 150 GeV. This choice of parameters leads to significant entropy 
production at temperatures immediately above the freeze-out surface, which is observed as a soft 
region in the speed of sound. 

Once the entropy density has been determined at all relevant temperatures, p{T) is evaluated via 
the integral 

rT 



p{T)=p{Tf)+ [ s{T')dT' 



(19) 



This procedure was employed to generate e, p, and s at a spacing of 1 MeV in temperature. The 
equation of state was then smoothed via a traditional cubic spline ^34] to ensure consistency and 
smoothness in the interpolation. 

The only remaining thermodynamic variables to be determined are related to the viscosity of the 
matter. It is generally thought that away from the cross-over temperature the ratio of shear viscosity 
to entropy density should be larger, although most models have fixed this ratio for all temperatures. 




T [GeV] 



FIG. 1. (color online) The squared speed of sound as a function of the temperature from the hadron resonance 
gas (red dotted), a fit to lattice data [3T] (blue dashed), and the merging of the two (black solid). The fit to 
the lattice data predicts entropy production at a lower temperature than observed in the hadron gas as seen in 
the lower speed of sound. Matching the entropy at T; then requires a soft region between Th and Tj. 
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FIG. 2. (color online) Snapshot of the energy density at r^s = and r = 5.52 fm/c along the x and y axes with 
and without Bjorkcn invariance for b = 5.70 fm. The central energy density falls more rapidly without Bjorken 
invariance, but the effect on the evolution diminishes in the periphery. 

However, viscous suppression of anisotropic flow at RHIC seems mostly to be driven by the value of 
the shear viscosity at the freeze-surface surface [35] although this may not be true at LHC energies 
|36| . In some cases, the low temperature regions of this model were found to be unstable with this 
choice and so at very low temperatures the shear viscosity was chosen to continuously vary with the 
energy density below some temperature well outside of the freeze-out surface (~ 120 MeV). The model 
was found to be insensitive to the choice of this temperature as long as the difference between it and 
the freeze-out temperature was more than 20 MeV. 

Finally, we take the relaxation time to be = This result can derived from linear response 
theory under the assumptions that the field loses correlation and that the entropy penalty for small 
viscous corrections is quadratic [T^ . 

VI. RESULTS 

We consider three centrality classes corresponding to three impact parameters. For the 0-5% 
centrality class we use b = 2.21 fm, b = 5.70 fm for 10-20%, and b = 7.23 for 20-30%, which are the 
values determined by the STAR collaboration using a Monte Carlo Glauber model [1]. The energy 
scale, Ko, fixes the maximum energy density in the most central collisions (b=2.21 fm) to be 19.5 
GeV/fm'^ at tq = 0.8 fm/c. The shear viscosity to entropy density ratio is taken to be 0.16, the initial 
value of the longitudinal part of the shear tensor vr^^ = —p/2, and the initial transverse flow parameter 
is taken as ^ = 0.5. 
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FIG. 3. (color online) The longitudinal velocity gradient at the origin equals l/r for a pure Bjorken expansion, 
so the ratio above is explicitly zero in the two-dimensional treatments. The velocity gradient varies from the 
Bjorken assumption in the three dimensional model on the order of 5%. 

We begin by demonstrating the effects of the finite longitudinal extent on the hydrodynamic evo- 
lution. Figure [2] shows the energy density at a fixed time along the transverse axes at mid-rapidity. 
The energy density near the center of the fireball falls somewhat faster in the full three-dimensional 
model. This was expected given that the additional longitudinal acceleration, which is neglected in a 
boost-invariant expansion, which leads to a slightly more rapid expansion and lower energy densities. 
However, the evolution in the periphery has changed surprisingly little. 

Figure [3] quantifies the extent to which the velocity gradient at the origin falls faster in the full 
three-dimensional model. Whereas the velocity gradient is dzU^ = 1/t in the boost-invariant picture, 
longitudinal acceleration increases the gradient. The increase from the Bjorken value of the longitu- 
dinal velocity shows little dependence on the impact parameter and is at the five percent level. This 
suggests that femtoscopic estimates of lifetime based on the boost-invariant assumption are a few 
percent low. This discrepancy should be somewhat larger for smaller viscosities, because an increase 
in transverse pressure should increase the longitudinal acceleration. 

Figure |3] suggests that the change in the energy density profile at rjs = is due to the presence 
non-trivial longitudinal expansion, and not due to changes in the transverse flow. In Figure |4j we 
turn our attention to this transverse flow under the same conditions as Figure [2} Whereas the energy 
density differed by up to 5% at midrapidity, we find that the collective velocity differs by a maximum 
of 1% for the same initial conditions. This implies that any significant modification to the elliptic 
flow at midrapidity will not come from changes in the dynamics at rjg = 0, though we cannot address 
changes to the initial conditions require to fit multiplicities or the effects of thermal smearing in this 

11 




FIG. 4. (color online) Snapshot of the collective at r^s = and r = 5.52 fm/c along the x and y axes with and 
without Bjorken invariance for b = 5.70 fm. The collective velocity at 77s = is found to be unaffected by the 
non-trivial longitudinal expansion. 

study. 

Figure [5] shows the difference in the location of the freeze-out surface along the x-axis at mid- 
rapidity. Due to the increased longitudinal velocity gradient in the three-dimensional calculation, the 
system breaks up slightly faster. All of these changes suggest only modest changes to the mid-rapidity 
observables predicted by boost-invariant viscous hydrodynamics. 



VII. CONCLUSION AND OUTLOOK 

We have presented a new simulation of (3-1-1) dimensional Israel-Stewart viscous hydrodynamics 
for relativistic heavy ion collisions geared toward predictions using average initial conditions. The two- 
dimensional version of this model was found to be highly consistent with two-dimensional calculations 
from other groups. The algorithm conserves T^^ at a high level, and in the limit of zero viscosity, 
entropy was conserved at better than the 0.1% level. 

Investigations of the effect of the non-trivial longitudinal expansion show the strongest effect at 
the level of 5%, as the system expands and disintegrates modestly faster than in corresponding two- 
dimensional calculations. The effects of the non-trivial expansion are stronger near the center of the 
fireball than in the periphery. Combined with very small changes to the transverse velocities, this 
causes the motion of the freeze-out surface to remain almost unchanged until later times when it falls 
more rapidly. This suggests that changes to important midrapidity observables like elliptic flow will 
not be larger than a few percent. 
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FIG. 5. (color online) The position of the freeze-out surface along the x-axis at mid-rapidity for three impact 
parameters. Incorporating the full three-dimensional motion reduces the lifetime of the source by about 5%. 
However, since much of the emission comes from the sudden, final collapse of the source the effect on emission 
may be modest. 



13 



These results suggest that one can ignore the effects of the full three-dimensional expansion if one 
is satisfied with predictions at the 5% level. Given that many approximations in the model are only 
trustworthy at the 5% level, these corrections can be considered modest, but cannot be classified as 
fully negligible. If one were to repeat the comparison at lower energies, the effects would certainly be 
greater, while at higher energies such as at the LHC, boost invariance should be significantly more 
justified. 

We intend to implement a freeze-out algorithm in the near future. Currently, the interface be- 
tween the cascade and the hydrodynamic model only functions in the two-dimensional version of the 
code. Once the three-dimensional interface is finished, comparison with data will be performed with 
observables at and away from mid-rapidity. 
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Appendix A: Evaluation of Local Derivatives 

For our set of integration variables, we choose a mix of quantities measured in the frame moving 
with the energy density and quantities measured in a frame undergoing a fixed expansion in spatial 
rapidity, which we refer to as the mesh frame. Beginning from the equations of motion as viewed in 
the matter frame, we must boost terms into the mesh frame. For instance, to track the velocity as 
observed in the mesh frame it is prudent to evaluate equations of motion in terms of spatial derivatives 
of the collective velocity in the mesh frame. In addition, one needs to write the equations of motion 
in terms of derivatives evaluated on a fixed grid of the mesh frame proper time. 

We introduce the general form of a boost from a frame that observes some velocity u^^ to one that 
observes any other velocity n^. In terms of these velocities and the metric tensor {g^^), the boost is: 

^ ' ^ l-t-tx-n ^ ' 

which fulfills the requirement that k^'^Uv = n'^- 

This boost formula can be used to evaluate local derivatives of local velocities. We take n to be 
the velocity as viewed in the matter frame, which is by definition n'^ = (1,0,0,0), and u to be the 
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velocity as viewed in the mesh frame. Likewise, we define the local derivative as 5^ = k^^di,. Then by 
one method, we can evaluate local derivatives of local velocities by simply evaluating: 



dan'' 



1 + 7 



On,u'^ — 5q7 



(A2) 
(A3) 



1 + 7 

which can be used to evaluate both spatial and temporal derivatives of the collective velocity. For 
instance, local spatial derivatives can be expressed 

,k 



di'nP = diw' + vJ'vP + 



u u 



1 + 7 



dkw> 



7(1 + 7) 



divp + u'vP + 



1 + 7 



(A4) 



where latin indices run over the spatial coordinates and repeated indices are summed over. The 
structure of our code is such that the coefficients associated with time derivatives of tabulated variables 
are separated from spatial derivatives, leading to the notation 



7(1 + 7) 



(A5) 



where the prime indicates the removal of those terms proportional to time derivatives of tabulated 
variables. 

Also among our variables are terms defined in the frame of the collective motion. The above 
construction cannot be used to calculate derivatives of these variables. Instead, we consider a small 
boost in the neighborhood of the boost from the mesh frame to the fluid frame. In our notation, this 
will be referred to as 6A.{u — 5u, n) where the choice of sign allows us to treat the differential element 
as naturally covariant (where the differential operator is normally naturally contravariant) . Reducing 
to only terms linear in the deviation of the fluid velocity: 



(5A^''(u - 5u, n) = g^"" + 2n'*(u'^ - Sv!") 



(u" - Sui" + n^'){u'' - Su'' + n'^) 



1 + 7 — (^7 



A(it, n) 



(A6) 



6j 



1+7 ' 1+7 (1 + 7)^ 

This method confirms the above results for local derivatives of the collective velocity by noting that 

+ 



1 + 7 

where we have frequently used the orthogonality relationship u^Sufj, = 0. 



(A7) 
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Now, since we wish to obtain corrections to derivatives of a tensor in the matter frame, a boost 
back to the frame of the matter is required. The boost to be appUed to the matter frame tensor is 
then: 



which can be appHed symmetrically to the local shear tensor: 



VT' 



M/9 



1+7 1+7 

where we have used the choice of frame occasionally (n^vr^'^ = 0). 

If both indices on the shear tensor are time-like, 5tt^^ = 0. For one time-like index and one space-like 
index (z), 

(57r°* = (5A°"Ac,^7r^^ + 5A^"Ac,/}7r^° 

= (5«a7r"^ - Y^^XaTT"^ = -<(5n°^ (AlO) 

where the final equivalence can be seen as a consistency check. Finally, for two spatial indices, 

5tt'^ = [du'UaTT'^^ + 5u^UaTT'^' - u'SUc.T^'^^ - SUc,!^"^'] (All) 

As an example, consider the simplified Israel-Stewart equation 

rM'^n + '^''' = '^[NS) (A12) 

After boosting the time derivative from the frame of the matter to the mesh frame via Bq = ^da + u^di, 
there remains a co- moving time derivative of a spatial component of the local shear tensor {'jdoTr'^-')- 
Simply integrating forward with respect to time would produce errors due to the changing reference 
frame in which tt is defined. Therefore, one amends this via 

doTT'^ + lu^vf'lT^^ + 'K^vf'u^ - il'uV^ - TT^'^u'^ilA (A13) 

1 + 7 L -I 

Now, our calculation is based on a fixed mesh in the expanding coordinate system (r, x, y, jy). Given 
this coordinate system, it only makes sense to choose our variables to be defined in the fixed frame of 
this coordinate system (the mesh frame). So when discussing the collective velocity u, we refer not to 
the collective velocity observed in the laboratory frame - the center of mass frame for the collision - 
but rather to the collective velocity observed by a comover with the longitudinal expansion. However, 
our equations of motion will be derived in the frame of the matter - a third, potentially unique Lorentz 
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frame. In light of this, our approach will again be to determine the small correction to be applied to 
derivatives in the mesh frame and boost this correction back into the matter frame. 

We begin by considering the velocity as observed in the mesh frame it^ = {-y , , u'^ , u^} . If we 
were to view this vector in the laboratory frame (wf), we would observe instead 

uf = {7 cosh 7/ + sinhr/, u^,u^ , coshry + 7sinhry} (^14) 

So, viewed from a frame at r/ = 5r?, 7 = 7 + u^6rj and = + ydr]. This means that if between mesh 
cells at slightly differing spatial rapidity one observes that «^(0) = u^{6r]), the difference in frame 
dictates that 9^u^(0) = 7. 

Now, any particular mesh frame may accurately view itself as being at 77 = 0, and its relationships 
to its neighbors at slightly differing spatial rapidity are preserved. So if one considered a small boost 
SA^'^ = ni^rf — ri^n'^ where n'* = {1,0,0,0} and ij^ = {0,0,0,— Srj}, such a boost would adjust the 
observed velocity as required: 

(5A>" = {n^'ria - jfno)u'^ = n^'iu ■ rf) - -fff = 6^'^u^Sr^ + S^'^'^Sri (A15) 

Boosting this into the frame of the matter, 

A^''{u,n)A''^iu,n)6Aap = A'^"A^^(nar7^ - VaUp) 

= (27n^ - «^)(77'^ + 2{u ■ r])n'' - + n'']) 

1 + 7 

-(r/^ + 2{u ■ 7])ni' - + rl^']){2^rf - u") (A16) 

1 + 7 

Applying this to the matter frame velocity, we obtain several results 

(AAM)°n" = 7(u ■ r?) - 7(u • r?) = (A17) 

(AA5A)>" = -u^'iu ■ rj) + ^^^''(^^ • v) = -^^u'Srj (A18) 



(AA(5A)^n" = -u^{u ■ rj) - j 



and a few results for the shear tensor 



1 + 7 



7 

1 + 7. 



Sri (A19) 



(AAM)^'*'^ = -u'^SriTr"' (A20) 
(AA(5A)>"^ = -{u'tt"' - Mavr"'^)(5ry = (n^vr^^ + w%J"^)5r? (A21) 

which will suffice as linear corrections to the shear tensor are always of these forms. 
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Appendix B: Equations of Motion 



Since the equations of motion here are for a set of variables not previously used, we include a 
sketch of the derivation of the equations of motion. In general, we will not include terms relating to 
the longitudinal expansion but will include time derivatives of the collective velocity related to the 
motion of the matter frame. 

In the frame of the matter, the conservation equations take the form 

doe + {e + P) din' + tt'^ diU^ = (Bl) 
(e + P) don' + tt'^ don^ + diP + djir'^ = (B2) 

Excluded here are any terms proportional to ra*, for instance n^diir'^ , which are explicitly zero in this 
frame. 

Separating out the time derivatives, but leaving alone {din^)' as defined in Appendix A for energy 
conservation: 

(e + PK . , ^^^^ _ {uW^)u^ ., 
7 7(1+7) 

+u'die + (e + P)diu' + Tr'^idin^y = 

where the second term of the first line comes from 7 and final two terms of the first line come from 
din^. 

Doing the same for momentum conservation: 

u'c^e + (e + P)^e - ^^^^Y^it^' + IT^'^u^ (B4) 



4 



k 



+ {'y — 1)'k'^u' + u^TT 



1 + 7 

+(e + P){don'y + TT'^idon^y + diP + ^^duP + {dji^'^y = 

1 + 7 



+ 7 

where several of the terms at the end of the first line come from djir'^ uWott'^ ■ 

In the frame of the matter, the Israel-Stewart equation becomes a relaxation equation: 

tJot^'^ + r,(4/3)7r^^5fcn^ = tt^^^^ - tt'^ (B5) 

which, for example, for the first projected shear stress element is: 

7^001 + u'diai + (4/3)aia^n'^ + ^ = ZR (^^^ _ 9^^^) (B6) 
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Extracting all of the time derivatives, we obtain 



2 .... 

7di - ^ ^ [(u^TT^-^ - n^vr^-') tt-' - u^Ti^'^if + u%-'^tt^] (B7) 



tA "^7(1 + 7)7 37 

where Sjai still contains corrections due to velocity gradients. 



Appendix C: Exponential Test 

One of the tests performed to validate the (3+1) dimensional code was the integration of a stationary 
exponential profile in energy density with length scale R: 



eo • e 



-x/R 



u = 



(CI) 



Since the energy gradients are then proportional to their values, the ideal hydrodynamic equations 
become linear and the collective velocity remains the same at every point. One may then solve for the 
transverse collective velocity in the presence of a Bjorken expansion, yielding a simple linear differential 
equation: 



u 



UQCI 



UqU 



+ 



1 



(C2) 



r R{1 + cj 

where the first term is from the Bjorken expansion and the second is from the exponential gradient. 
In the case where there is no Bjorken Expansion, the first term is zero and the term equation has an 
analytical solution: 

R^llc2/ = ^ + sinh-n^) - (c's - l)nx/n^] (C3) 

In the presence of the Bjorken expansion the equation is not separable and no analytical solution 
was pursued. The same computation for a longitudinal exponential profile yields the same differential 
equation with the assignment R — t- ct^t in the solution. The solutions for the these three cases are 
presented in Fig. [6j 

One can then solve for the evolution of the energy density, which can be written in terms of the 
velocity as 



Uq 



^2(1 



+ 1 



(1 



[l + ci)uo 



R 



(C4) 



The general case can be easily integrated numerically, but there is a semi-analytic solution for the 
case with no Bjorken expansion. One exploites a change of variables to arrive at an expression for the 



19 




FIG. 6. (color online) The collective velocity at all points for a exponential energy density profile with no 
Bjorken expansion (blue squares), a transverse profile with Bjorken expansion (red stars), and a longitudinal 
profile with Bjorken expansion (green triangles). In each case the velocity calculated by the full code agrees 
with the correct value. 
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energy density given the velocity 



du- 



dt de 
du dt 
1 - 



eo • exp 



2c2 



(C5) 



The results for a system with no expansion, a transverse system with expansion, and a longitudinal 
system with expansion are shown in Fig. [7j 
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